We start by loading in our data tables

require(mgcv)
bologna <- read.csv("bologna_full.csv")
bologna <- na.omit(bologna)

modena <- read.csv("modena_full.csv")
modena <- na.omit(modena)

florence <- read.csv("florence_full.csv")
florence <- na.omit(florence)

rome <- read.csv("rome_full.csv")
rome <- na.omit(rome)

num_iter <- 30

First lets look at effects of NPI

Bologna:

for(i in 0:num_iter) {
num_lag <- i
bologna <- read.csv("bologna_full.csv")
bologna$Daily.1 <- c(bologna$Daily[(1+num_lag):(nrow(bologna))], rep(NA, num_lag))
bologna <- na.omit(bologna)

infections.gam <- gam(Daily.1 ~ s(NPI, k=5), family=poisson(link="log"), data=bologna)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
summary(infections.gam)
print(cat("Day:", i, "R-sq", summary(infections.gam)$r.sq))

}

Modena:

for(i in 0:num_iter) {
num_lag <- i
modena <- read.csv("modena_full.csv")
modena$Daily.1 <- c(modena$Daily[(1+num_lag):(nrow(modena))], rep(NA, num_lag))
modena <- na.omit(modena)

infections.gam <- gam(Daily.1 ~  s(NPI, k=5), family=poisson(link="log"), data=modena)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
summary(infections.gam)
print(cat("Day:", i, "R-sq", summary(infections.gam)$r.sq))

}

Florence:

for(i in 0:num_iter) {
num_lag <- i
florence <- read.csv("florence_full.csv")
florence$Daily.1 <- c(florence$Daily[(1+num_lag):(nrow(florence))], rep(NA, num_lag))
florence <- na.omit(florence)

infections.gam <- gam(Daily.1 ~ s(NPI, k=5), family=poisson(link="log"), data=florence)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
summary(infections.gam)
print(cat("Day:", i, "R-sq", summary(infections.gam)$r.sq))

}

Rome:

for(i in 0:num_iter) {
num_lag <- i
rome <- read.csv("rome_full.csv")
rome$Daily.1 <- c(rome$Daily[(num_lag+1):(nrow(rome))], rep(NA, num_lag))
rome <- na.omit(rome)

infections.gam <- gam(Daily.1 ~  s(NPI, k=5), family=poisson(link="log"), data=rome)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
print(cat("Day:", i, "R-sq", summary(infections.gam)$r.sq))
}

Now lets look at Temp/Humidity

Bologna:

for(i in 0:num_iter) {
num_lag <- i
bologna <- read.csv("bologna_full.csv")
bologna$Daily.1 <- c(bologna$Daily[(1+num_lag):(nrow(bologna))], rep(NA, num_lag))
bologna <- na.omit(bologna)

infections.gam <- gam(Daily.1 ~ s(Temperature) + s(Humidity), family=poisson(link="log"), data=bologna)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
summary(infections.gam)
print(cat("Day:", i, "R-sq", summary(infections.gam)$r.sq))

}
Day: 0 R-sq 0.4896344NULL

Day: 1 R-sq 0.4633273NULL

Day: 2 R-sq 0.4918461NULL

Day: 3 R-sq 0.5703598NULL

Day: 4 R-sq 0.5318496NULL

Day: 5 R-sq 0.4716657NULL

Day: 6 R-sq 0.4663372NULL

Day: 7 R-sq 0.5765711NULL

Day: 8 R-sq 0.5500971NULL

Day: 9 R-sq 0.5945654NULL

Day: 10 R-sq 0.6034324NULL

Day: 11 R-sq 0.591467NULL

Day: 12 R-sq 0.6489848NULL

Day: 13 R-sq 0.6266869NULL

Day: 14 R-sq 0.6410671NULL

Day: 15 R-sq 0.6286455NULL

Day: 16 R-sq 0.6710816NULL

Day: 17 R-sq 0.6710324NULL

Day: 18 R-sq 0.6733654NULL

Day: 19 R-sq 0.7549484NULL

Day: 20 R-sq 0.7704205NULL

Day: 21 R-sq 0.7335032NULL

Day: 22 R-sq 0.6961922NULL

Day: 23 R-sq 0.7280413NULL

Day: 24 R-sq 0.7578183NULL

Day: 25 R-sq 0.7361595NULL

Day: 26 R-sq 0.6891041NULL

Day: 27 R-sq 0.6511969NULL

Day: 28 R-sq 0.6365822NULL

Day: 29 R-sq 0.617207NULL

Day: 30 R-sq 0.5663472NULL

Modena:

for(i in 0:num_iter) {
num_lag <- i
modena <- read.csv("modena_full.csv")
modena$Daily.1 <- c(modena$Daily[(1+num_lag):(nrow(modena))], rep(NA, num_lag))
modena <- na.omit(modena)

infections.gam <- gam(Daily.1 ~  s(Temperature) + s(Humidity), family=poisson(link="log"), data=modena)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
summary(infections.gam)
print(cat("Day:", i, "R-sq", summary(infections.gam)$r.sq))

}
Day: 0 R-sq 0.5380692NULL

Day: 1 R-sq 0.4603939NULL

Day: 2 R-sq 0.5502608NULL

Day: 3 R-sq 0.5257741NULL

Day: 4 R-sq 0.495446NULL

Day: 5 R-sq 0.412964NULL

Day: 6 R-sq 0.4962353NULL

Day: 7 R-sq 0.553577NULL

Day: 8 R-sq 0.5091717NULL

Day: 9 R-sq 0.4840658NULL

Day: 10 R-sq 0.6148442NULL

Day: 11 R-sq 0.6935744NULL

Day: 12 R-sq 0.601654NULL

Day: 13 R-sq 0.6328875NULL

Day: 14 R-sq 0.6602823NULL

Day: 15 R-sq 0.6418581NULL

Day: 16 R-sq 0.7038114NULL

Day: 17 R-sq 0.6879685NULL

Day: 18 R-sq 0.6928742NULL

Day: 19 R-sq 0.7349044NULL

Day: 20 R-sq 0.6730137NULL

Day: 21 R-sq 0.6289489NULL

Day: 22 R-sq 0.6293349NULL

Day: 23 R-sq 0.5951238NULL

Day: 24 R-sq 0.6252658NULL

Day: 25 R-sq 0.691863NULL

Day: 26 R-sq 0.5791828NULL

Day: 27 R-sq 0.5838233NULL

Day: 28 R-sq 0.5648786NULL

Day: 29 R-sq 0.6451174NULL

Day: 30 R-sq 0.6263023NULL

Florence:

for(i in 0:num_iter) {
num_lag <- i
florence <- read.csv("florence_full.csv")
florence$Daily.1 <- c(florence$Daily[(1+num_lag):(nrow(florence))], rep(NA, num_lag))
florence <- na.omit(florence)

infections.gam <- gam(Daily.1 ~ s(Temperature) + s(Humidity), family=poisson(link="log"), data=florence)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
summary(infections.gam)
print(cat("Day:", i, "R-sq", summary(infections.gam)$r.sq))

}
Day: 0 R-sq 0.4728276NULL

Day: 1 R-sq 0.3744223NULL

Day: 2 R-sq 0.4542527NULL

Day: 3 R-sq 0.5526952NULL

Day: 4 R-sq 0.4970105NULL

Day: 5 R-sq 0.4053944NULL

Day: 6 R-sq 0.3662819NULL

Day: 7 R-sq 0.6391539NULL

Day: 8 R-sq 0.664867NULL

Day: 9 R-sq 0.6566028NULL

Day: 10 R-sq 0.5371588NULL

Day: 11 R-sq 0.5609689NULL

Day: 12 R-sq 0.5139698NULL

Day: 13 R-sq 0.4925803NULL

Day: 14 R-sq 0.4693117NULL

Day: 15 R-sq 0.5636821NULL

Day: 16 R-sq 0.6640505NULL

Day: 17 R-sq 0.6977326NULL

Day: 18 R-sq 0.5019629NULL

Day: 19 R-sq 0.5811804NULL

Day: 20 R-sq 0.6194071NULL

Day: 21 R-sq 0.5790887NULL

Day: 22 R-sq 0.5453216NULL

Day: 23 R-sq 0.6266673NULL

Day: 24 R-sq 0.5540829NULL

Day: 25 R-sq 0.6299241NULL

Day: 26 R-sq 0.6185027NULL

Day: 27 R-sq 0.5877647NULL

Day: 28 R-sq 0.5080166NULL

Day: 29 R-sq 0.6840771NULL

Day: 30 R-sq 0.7033451NULL

Rome:

for(i in 0:num_iter) {
num_lag <- i
rome <- read.csv("rome_full.csv")
rome$Daily.1 <- c(rome$Daily[(num_lag+1):(nrow(rome))], rep(NA, num_lag))
rome <- na.omit(rome)

infections.gam <- gam(Daily.1 ~  s(Temperature) + s(Humidity), family=poisson(link="log"), data=rome)
plot(infections.gam,scale=0,se=2, shade=TRUE,pages=1)
title(main = i)
print(cat("Day:", i, "R-sq", summary(infections.gam)$r.sq))
}
Day: 0 R-sq 0.4202767NULL

Day: 1 R-sq 0.4343473NULL

Day: 2 R-sq 0.4526714NULL

Day: 3 R-sq 0.4636081NULL

Day: 4 R-sq 0.4757067NULL

Day: 5 R-sq 0.5211681NULL

Day: 6 R-sq 0.5508436NULL

Day: 7 R-sq 0.5854002NULL

Day: 8 R-sq 0.6340703NULL

Day: 9 R-sq 0.6287562NULL

Day: 10 R-sq 0.5980017NULL

Day: 11 R-sq 0.6214937NULL

Day: 12 R-sq 0.6830744NULL

Day: 13 R-sq 0.7555303NULL

Day: 14 R-sq 0.7050403NULL

Day: 15 R-sq 0.7038004NULL

Day: 16 R-sq 0.7667924NULL

Day: 17 R-sq 0.806356NULL

Day: 18 R-sq 0.7858903NULL

Day: 19 R-sq 0.7727713NULL

Day: 20 R-sq 0.7579232NULL

Day: 21 R-sq 0.774631NULL

Day: 22 R-sq 0.7890732NULL

Day: 23 R-sq 0.7949118NULL

Day: 24 R-sq 0.7691577NULL

Day: 25 R-sq 0.7223769NULL

Day: 26 R-sq 0.7041938NULL

Day: 27 R-sq 0.7103775NULL

Day: 28 R-sq 0.7876148NULL

Day: 29 R-sq 0.7265889NULL

Day: 30 R-sq 0.7234478NULL

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKV2Ugc3RhcnQgYnkgbG9hZGluZyBpbiBvdXIgZGF0YSB0YWJsZXMKCmBgYHtyfQpyZXF1aXJlKG1nY3YpCmJvbG9nbmEgPC0gcmVhZC5jc3YoImJvbG9nbmFfZnVsbC5jc3YiKQpib2xvZ25hIDwtIG5hLm9taXQoYm9sb2duYSkKCm1vZGVuYSA8LSByZWFkLmNzdigibW9kZW5hX2Z1bGwuY3N2IikKbW9kZW5hIDwtIG5hLm9taXQobW9kZW5hKQoKZmxvcmVuY2UgPC0gcmVhZC5jc3YoImZsb3JlbmNlX2Z1bGwuY3N2IikKZmxvcmVuY2UgPC0gbmEub21pdChmbG9yZW5jZSkKCnJvbWUgPC0gcmVhZC5jc3YoInJvbWVfZnVsbC5jc3YiKQpyb21lIDwtIG5hLm9taXQocm9tZSkKCm51bV9pdGVyIDwtIDMwCmBgYAoKRmlyc3QgbGV0cyBsb29rIGF0IGVmZmVjdHMgb2YgTlBJCgoKQm9sb2duYToKYGBge3J9CmZvcihpIGluIDA6bnVtX2l0ZXIpIHsKbnVtX2xhZyA8LSBpCmJvbG9nbmEgPC0gcmVhZC5jc3YoImJvbG9nbmFfZnVsbC5jc3YiKQpib2xvZ25hJERhaWx5LjEgPC0gYyhib2xvZ25hJERhaWx5WygxK251bV9sYWcpOihucm93KGJvbG9nbmEpKV0sIHJlcChOQSwgbnVtX2xhZykpCmJvbG9nbmEgPC0gbmEub21pdChib2xvZ25hKQoKaW5mZWN0aW9ucy5nYW0gPC0gZ2FtKERhaWx5LjEgfiBzKE5QSSwgaz01KSwgZmFtaWx5PXBvaXNzb24obGluaz0ibG9nIiksIGRhdGE9Ym9sb2duYSkKcGxvdChpbmZlY3Rpb25zLmdhbSxzY2FsZT0wLHNlPTIsIHNoYWRlPVRSVUUscGFnZXM9MSkKdGl0bGUobWFpbiA9IGkpCnN1bW1hcnkoaW5mZWN0aW9ucy5nYW0pCnByaW50KGNhdCgiRGF5OiIsIGksICJSLXNxIiwgc3VtbWFyeShpbmZlY3Rpb25zLmdhbSkkci5zcSkpCgp9CmBgYAoKTW9kZW5hOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKbW9kZW5hIDwtIHJlYWQuY3N2KCJtb2RlbmFfZnVsbC5jc3YiKQptb2RlbmEkRGFpbHkuMSA8LSBjKG1vZGVuYSREYWlseVsoMStudW1fbGFnKToobnJvdyhtb2RlbmEpKV0sIHJlcChOQSwgbnVtX2xhZykpCm1vZGVuYSA8LSBuYS5vbWl0KG1vZGVuYSkKCmluZmVjdGlvbnMuZ2FtIDwtIGdhbShEYWlseS4xIH4gIHMoTlBJLCBrPTUpLCBmYW1pbHk9cG9pc3NvbihsaW5rPSJsb2ciKSwgZGF0YT1tb2RlbmEpCnBsb3QoaW5mZWN0aW9ucy5nYW0sc2NhbGU9MCxzZT0yLCBzaGFkZT1UUlVFLHBhZ2VzPTEpCnRpdGxlKG1haW4gPSBpKQpzdW1tYXJ5KGluZmVjdGlvbnMuZ2FtKQpwcmludChjYXQoIkRheToiLCBpLCAiUi1zcSIsIHN1bW1hcnkoaW5mZWN0aW9ucy5nYW0pJHIuc3EpKQoKfQpgYGAKCkZsb3JlbmNlOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKZmxvcmVuY2UgPC0gcmVhZC5jc3YoImZsb3JlbmNlX2Z1bGwuY3N2IikKZmxvcmVuY2UkRGFpbHkuMSA8LSBjKGZsb3JlbmNlJERhaWx5WygxK251bV9sYWcpOihucm93KGZsb3JlbmNlKSldLCByZXAoTkEsIG51bV9sYWcpKQpmbG9yZW5jZSA8LSBuYS5vbWl0KGZsb3JlbmNlKQoKaW5mZWN0aW9ucy5nYW0gPC0gZ2FtKERhaWx5LjEgfiBzKE5QSSwgaz01KSwgZmFtaWx5PXBvaXNzb24obGluaz0ibG9nIiksIGRhdGE9ZmxvcmVuY2UpCnBsb3QoaW5mZWN0aW9ucy5nYW0sc2NhbGU9MCxzZT0yLCBzaGFkZT1UUlVFLHBhZ2VzPTEpCnRpdGxlKG1haW4gPSBpKQpzdW1tYXJ5KGluZmVjdGlvbnMuZ2FtKQpwcmludChjYXQoIkRheToiLCBpLCAiUi1zcSIsIHN1bW1hcnkoaW5mZWN0aW9ucy5nYW0pJHIuc3EpKQoKfQpgYGAKCgpSb21lOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKcm9tZSA8LSByZWFkLmNzdigicm9tZV9mdWxsLmNzdiIpCnJvbWUkRGFpbHkuMSA8LSBjKHJvbWUkRGFpbHlbKG51bV9sYWcrMSk6KG5yb3cocm9tZSkpXSwgcmVwKE5BLCBudW1fbGFnKSkKcm9tZSA8LSBuYS5vbWl0KHJvbWUpCgppbmZlY3Rpb25zLmdhbSA8LSBnYW0oRGFpbHkuMSB+ICBzKE5QSSwgaz01KSwgZmFtaWx5PXBvaXNzb24obGluaz0ibG9nIiksIGRhdGE9cm9tZSkKcGxvdChpbmZlY3Rpb25zLmdhbSxzY2FsZT0wLHNlPTIsIHNoYWRlPVRSVUUscGFnZXM9MSkKdGl0bGUobWFpbiA9IGkpCnByaW50KGNhdCgiRGF5OiIsIGksICJSLXNxIiwgc3VtbWFyeShpbmZlY3Rpb25zLmdhbSkkci5zcSkpCn0KYGBgCgoKCgoKCgoKCgoKCgpOb3cgbGV0cyBsb29rIGF0IFRlbXAvSHVtaWRpdHkKCgpCb2xvZ25hOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKYm9sb2duYSA8LSByZWFkLmNzdigiYm9sb2duYV9mdWxsLmNzdiIpCmJvbG9nbmEkRGFpbHkuMSA8LSBjKGJvbG9nbmEkRGFpbHlbKDErbnVtX2xhZyk6KG5yb3coYm9sb2duYSkpXSwgcmVwKE5BLCBudW1fbGFnKSkKYm9sb2duYSA8LSBuYS5vbWl0KGJvbG9nbmEpCgppbmZlY3Rpb25zLmdhbSA8LSBnYW0oRGFpbHkuMSB+IHMoVGVtcGVyYXR1cmUpICsgcyhIdW1pZGl0eSksIGZhbWlseT1wb2lzc29uKGxpbms9ImxvZyIpLCBkYXRhPWJvbG9nbmEpCnBsb3QoaW5mZWN0aW9ucy5nYW0sc2NhbGU9MCxzZT0yLCBzaGFkZT1UUlVFLHBhZ2VzPTEpCnRpdGxlKG1haW4gPSBpKQpzdW1tYXJ5KGluZmVjdGlvbnMuZ2FtKQpwcmludChjYXQoIkRheToiLCBpLCAiUi1zcSIsIHN1bW1hcnkoaW5mZWN0aW9ucy5nYW0pJHIuc3EpKQoKfQpgYGAKCk1vZGVuYToKYGBge3J9CmZvcihpIGluIDA6bnVtX2l0ZXIpIHsKbnVtX2xhZyA8LSBpCm1vZGVuYSA8LSByZWFkLmNzdigibW9kZW5hX2Z1bGwuY3N2IikKbW9kZW5hJERhaWx5LjEgPC0gYyhtb2RlbmEkRGFpbHlbKDErbnVtX2xhZyk6KG5yb3cobW9kZW5hKSldLCByZXAoTkEsIG51bV9sYWcpKQptb2RlbmEgPC0gbmEub21pdChtb2RlbmEpCgppbmZlY3Rpb25zLmdhbSA8LSBnYW0oRGFpbHkuMSB+ICBzKFRlbXBlcmF0dXJlKSArIHMoSHVtaWRpdHkpLCBmYW1pbHk9cG9pc3NvbihsaW5rPSJsb2ciKSwgZGF0YT1tb2RlbmEpCnBsb3QoaW5mZWN0aW9ucy5nYW0sc2NhbGU9MCxzZT0yLCBzaGFkZT1UUlVFLHBhZ2VzPTEpCnRpdGxlKG1haW4gPSBpKQpzdW1tYXJ5KGluZmVjdGlvbnMuZ2FtKQpwcmludChjYXQoIkRheToiLCBpLCAiUi1zcSIsIHN1bW1hcnkoaW5mZWN0aW9ucy5nYW0pJHIuc3EpKQoKfQpgYGAKCkZsb3JlbmNlOgpgYGB7cn0KZm9yKGkgaW4gMDpudW1faXRlcikgewpudW1fbGFnIDwtIGkKZmxvcmVuY2UgPC0gcmVhZC5jc3YoImZsb3JlbmNlX2Z1bGwuY3N2IikKZmxvcmVuY2UkRGFpbHkuMSA8LSBjKGZsb3JlbmNlJERhaWx5WygxK251bV9sYWcpOihucm93KGZsb3JlbmNlKSldLCByZXAoTkEsIG51bV9sYWcpKQpmbG9yZW5jZSA8LSBuYS5vbWl0KGZsb3JlbmNlKQoKaW5mZWN0aW9ucy5nYW0gPC0gZ2FtKERhaWx5LjEgfiBzKFRlbXBlcmF0dXJlKSArIHMoSHVtaWRpdHkpLCBmYW1pbHk9cG9pc3NvbihsaW5rPSJsb2ciKSwgZGF0YT1mbG9yZW5jZSkKcGxvdChpbmZlY3Rpb25zLmdhbSxzY2FsZT0wLHNlPTIsIHNoYWRlPVRSVUUscGFnZXM9MSkKdGl0bGUobWFpbiA9IGkpCnN1bW1hcnkoaW5mZWN0aW9ucy5nYW0pCnByaW50KGNhdCgiRGF5OiIsIGksICJSLXNxIiwgc3VtbWFyeShpbmZlY3Rpb25zLmdhbSkkci5zcSkpCgp9CmBgYAoKClJvbWU6CmBgYHtyfQpmb3IoaSBpbiAwOm51bV9pdGVyKSB7Cm51bV9sYWcgPC0gaQpyb21lIDwtIHJlYWQuY3N2KCJyb21lX2Z1bGwuY3N2IikKcm9tZSREYWlseS4xIDwtIGMocm9tZSREYWlseVsobnVtX2xhZysxKToobnJvdyhyb21lKSldLCByZXAoTkEsIG51bV9sYWcpKQpyb21lIDwtIG5hLm9taXQocm9tZSkKCmluZmVjdGlvbnMuZ2FtIDwtIGdhbShEYWlseS4xIH4gIHMoVGVtcGVyYXR1cmUpICsgcyhIdW1pZGl0eSksIGZhbWlseT1wb2lzc29uKGxpbms9ImxvZyIpLCBkYXRhPXJvbWUpCnBsb3QoaW5mZWN0aW9ucy5nYW0sc2NhbGU9MCxzZT0yLCBzaGFkZT1UUlVFLHBhZ2VzPTEpCnRpdGxlKG1haW4gPSBpKQpwcmludChjYXQoIkRheToiLCBpLCAiUi1zcSIsIHN1bW1hcnkoaW5mZWN0aW9ucy5nYW0pJHIuc3EpKQp9CmBgYAo=